home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / iritsm3s.zip / PRIMITIV.C < prev    next >
C/C++ Source or Header  |  1992-01-25  |  52KB  |  1,327 lines
  1. /*****************************************************************************
  2. *   "Irit" - the 3d polygonal solid modeller.                     *
  3. *                                         *
  4. * Written by:  Gershon Elber                Ver 0.2, Mar. 1990   *
  5. ******************************************************************************
  6. *   Module to generate the geometric primitives defined in the system. The   *
  7. * primitives currently defined are:                         *
  8. * 1. BOX - main planes parallel box.                         *
  9. * 2. GBOX - generalized box - 6 arbitrary planes.                 *
  10. * 3. CYLIN - cylinder with any main direction.                     *
  11. * 4. CONE, CONE2 - cone with any main direction (two bases).             *
  12. * 5. SPHERE                                     *
  13. * 6. TORUS - with any main direction.                         *
  14. * 7. PLANE - non closed, single polygon object: circle with resolution edges *
  15. * 8. POLY - directly define single polygon object by specifing its vertices. *
  16. *   In addition, the following lower level operations are defined to create  *
  17. * objects - EXTRUDE, and SURFREV, both require a polygon and a vector to     *
  18. * extrude/rotate the polygon along.                         *
  19. *****************************************************************************/
  20.  
  21. #include <stdio.h>
  22. #include <math.h>
  23. #include "program.h"
  24. #include "allocate.h"
  25. #include "attribut.h"
  26. #include "convex.h"
  27. #include "geomat3d.h"
  28. #include "graphgen.h"
  29. #include "objects.h"
  30. #include "primitiv.h"
  31. #include "windows.h"
  32.  
  33. #define    MIN_RESOLUTION 4
  34.  
  35. static PolygonStruct *GenInsidePoly(PolygonStruct *Pl);
  36. static PolygonStruct *GenPolygon4Vrtx(VectorType V1, VectorType V2,
  37.     VectorType V3, VectorType V4, VectorType Vin, PolygonStruct *Pnext);
  38. static PolygonStruct *GenPolygon3Vrtx(VectorType V1, VectorType V2,
  39.             VectorType V3, VectorType Vin, PolygonStruct *Pnext);
  40. static void GenTransformMatrix(MatrixType Mat, VectorType Trans,
  41.                         VectorType Dir, RealType Scale);
  42. static void UpdateVertexNormal(NormalType Normal, PointType Pt, PointType InPt,
  43.                     int Perpendicular, PointType PerpPt);
  44.  
  45. /*****************************************************************************
  46. *   Routine to create a BOX geometric object defined by Pt - the minimun     *
  47. * 3d point, and Width - Dx Dy & Dz vector.        4             *
  48. * Order of vertices is as                           5       7             *
  49. * follows in the picture:                           |   6   |             *
  50. *                            |   |   |             *
  51. * (Note vertex 0 is hidden behind edge 2-6)        |    |   |             *
  52. *                            1   |   3                *
  53. *                            2             *
  54. *****************************************************************************/
  55. ObjectStruct * GenBOXObject(VectorType Pt, RealType *WidthX,
  56.                     RealType *WidthY, RealType *WidthZ)
  57. {
  58.     VectorType Dir1, Dir2, Dir3;
  59.  
  60.     PT_CLEAR(Dir1);    Dir1[0] = (*WidthX);   /* Prepare direction vectors. */
  61.     PT_CLEAR(Dir2);    Dir2[1] = (*WidthY);       /* Parallel to main axes. */
  62.     PT_CLEAR(Dir3);    Dir3[2] = (*WidthZ);           /* For GBOX call. */
  63.  
  64.     return GenGBOXObject(Pt, Dir1, Dir2, Dir3);
  65. }
  66.  
  67. /*****************************************************************************
  68. *   Routine to create a GBOX geometric object defined by Pt - the minimun    *
  69. * 3d point, and 3 direction Vectors Dir1, Dir2, Dir3. If two of the         *
  70. * direction vectors are parallel the GBOX converges to zero volume. A NULL   *
  71. * pointer is returned in that case!                         *
  72. *                             4             *
  73. * Order of vertices is as                           5       7             *
  74. * follows in the picture:                           |   6   |             *
  75. *                            |   |   |             *
  76. * (Note vertex 0 is hidden behind edge 2-6)        |    |   |             *
  77. *                            1   |   3                *
  78. *                            2             *
  79. *****************************************************************************/
  80. ObjectStruct * GenGBOXObject(VectorType Pt,
  81.             VectorType Dir1, VectorType Dir2, VectorType Dir3)
  82. {
  83.     int i;
  84.     VectorType Temp;
  85.     VectorType V[8];                  /* Hold 8 vertices of BOX. */
  86.     VertexStruct *PVertex;
  87.     PolygonStruct *PPolygon;
  88.     ObjectStruct *PBox;
  89.  
  90.     VecCrossProd(Temp, Dir1, Dir2);
  91.     if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
  92.     VecCrossProd(Temp, Dir2, Dir3);
  93.     if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
  94.     VecCrossProd(Temp, Dir3, Dir1);
  95.     if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
  96.  
  97.     /* Also the 0..7 sequence is binary decoded such that bit 0 is Dir1, */
  98.     /* bit 1 Dir2, and bit 2 is Dir3 increment:                 */
  99.     for (i = 0; i < 8; i++) {
  100.     PT_COPY(V[i], Pt);
  101.  
  102.     if (i & 1) { PT_ADD(V[i], V[i], Dir1); }
  103.     if (i & 2) { PT_ADD(V[i], V[i], Dir2); }
  104.     if (i & 4) { PT_ADD(V[i], V[i], Dir3); }
  105.     }
  106.  
  107.     PBox = GenPolyObject("", NULL, NULL); /* Generate the BOX object itself: */
  108.  
  109.     /* And generate the 6 polygons (Bottom, top and 4 sides in this order):  */
  110.     PBox -> U.Pl.P = GenPolygon4Vrtx(V[0], V[1], V[3], V[2], V[4], PBox -> U.Pl.P);
  111.     PBox -> U.Pl.P = GenPolygon4Vrtx(V[6], V[7], V[5], V[4], V[0], PBox -> U.Pl.P);
  112.     PBox -> U.Pl.P = GenPolygon4Vrtx(V[4], V[5], V[1], V[0], V[2], PBox -> U.Pl.P);
  113.     PBox -> U.Pl.P = GenPolygon4Vrtx(V[5], V[7], V[3], V[1], V[0], PBox -> U.Pl.P);
  114.     PBox -> U.Pl.P = GenPolygon4Vrtx(V[7], V[6], V[2], V[3], V[1], PBox -> U.Pl.P);
  115.     PBox -> U.Pl.P = GenPolygon4Vrtx(V[6], V[4], V[0], V[2], V[3], PBox -> U.Pl.P);
  116.  
  117.     /* Update the vertices normals using the polygon plane equation: */
  118.     for (PPolygon = PBox -> U.Pl.P;
  119.      PPolygon != NULL;
  120.      PPolygon = PPolygon -> Pnext) {
  121.     PVertex = PPolygon -> V;
  122.     do {
  123.         PT_COPY(PVertex -> Normal, PPolygon -> Plane);
  124.         PVertex = PVertex -> Pnext;
  125.     }
  126.     while (PVertex != PPolygon -> V);
  127.     }
  128.  
  129.  
  130.     SetObjectColor(PBox, GlblPrimColor);       /* Set its default color. */
  131.  
  132.     return PBox;
  133. }
  134.  
  135. /*****************************************************************************
  136. *  Routine to fetch the resolution parameter from the RESOLUTION object.     *
  137. *  If ClipToMin TRUE, the Resolution is clipped to not be below MIN_RES.     *
  138. *****************************************************************************/
  139. int GetResolution(int ClipToMin)
  140. {
  141.     int Resolution;
  142.     ObjectStruct *PObj = GetObject("RESOLUTION");
  143.  
  144.     if (PObj == NULL || !IS_NUM_OBJ(PObj)) {
  145.     WndwInputWindowPutStr("No numeric object name RESOLUTION is defined");
  146.     Resolution = DEFAULT_RESOLUTION;
  147.     }
  148.     else
  149.     Resolution = ClipToMin ? MAX(((int) (PObj -> U.R)), MIN_RESOLUTION)
  150.                    : (int) (PObj -> U.R);
  151.  
  152.     Resolution = (Resolution / 2) * 2;        /* Make sure its an even number. */
  153.  
  154.     return Resolution;
  155. }
  156.  
  157. /*****************************************************************************
  158. *   Routine to create a CONE geometric object defined by Pt - the base       *
  159. * 3d center point, Dir - the cone direction and length, and base radius R.   *
  160. *****************************************************************************/
  161. ObjectStruct *GenCONEObject(VectorType Pt, VectorType Dir, RealType *R)
  162. {
  163.     int i, Resolution;
  164.     RealType Angle, AngleStep;
  165.     PointType LastCirclePt, CirclePt, ApexPt;
  166.     NormalType LastCircleNrml, CircleNrml, ApexNrml;
  167.     MatrixType Mat;
  168.     VertexStruct *VBase, *PVertex;
  169.     PolygonStruct *PBase;
  170.     ObjectStruct *PCone;
  171.  
  172.     Resolution = GetResolution(TRUE);     /* Get refinement factor of object. */
  173.  
  174.     GenTransformMatrix(Mat, Pt, Dir, *R);     /* Transform from unit circle. */
  175.  
  176.     PT_COPY(ApexPt, Pt);           /* Find the apex point: Pt + Dir. */
  177.     PT_ADD(ApexPt, ApexPt, Dir);
  178.     PT_NORMALIZE(Dir);
  179.  
  180.     PCone = GenPolyObject("", NULL, NULL);   /* Gen. the CONE object itself: */
  181.     /* Also allocate the base polygon header with first vertex on it: */
  182.     PBase = AllocPolygon(0, 0, VBase = AllocVertex(0, 0, NULL, NULL), NULL);
  183.  
  184.     LastCirclePt[0] = 1.0;        /* First point is allways Angle = 0. */
  185.     LastCirclePt[1] = 0.0;
  186.     LastCirclePt[2] = 0.0;
  187.     MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);
  188.  
  189.     UpdateVertexNormal(LastCircleNrml, LastCirclePt, Pt, TRUE, ApexPt);
  190.  
  191.     PT_COPY(VBase -> Pt, LastCirclePt);  /* Update first pt in base polygon. */
  192.     PT_COPY(VBase -> Normal, Dir);
  193.  
  194.     AngleStep = M_PI * 2 / Resolution;
  195.  
  196.     for (i = 1; i <= Resolution; i++) {          /* Pass the whole base circle. */
  197.     Angle = AngleStep * i;             /* Prevent from additive error. */
  198.  
  199.     CirclePt[0] = cos(Angle);
  200.     CirclePt[1] = sin(Angle);
  201.     CirclePt[2] = 0.0;
  202.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  203.  
  204.      UpdateVertexNormal(CircleNrml, CirclePt, Pt, TRUE, ApexPt);
  205.  
  206.     PCone -> U.Pl.P = GenPolygon3Vrtx(LastCirclePt, ApexPt,
  207.                           CirclePt, Pt, PCone -> U.Pl.P);
  208.  
  209.     /* Update the normals for this cone side polygon vertices: */
  210.     PVertex = PCone -> U.Pl.P -> V;
  211.     PT_COPY(PVertex -> Normal, LastCircleNrml);
  212.     PVertex = PVertex -> Pnext;
  213.     /* The apex normal is the average of the two base vertices: */
  214.     PT_ADD(ApexNrml, CircleNrml, LastCircleNrml);
  215.     PT_NORMALIZE(ApexNrml);
  216.     PT_COPY(PVertex -> Normal, ApexNrml);
  217.     PVertex = PVertex -> Pnext;
  218.     PT_COPY(PVertex -> Normal, CircleNrml);
  219.  
  220.     /* And add this vertex to base polygon: */
  221.     if (i == Resolution)           /* Its last point - make it circular. */
  222.         VBase -> Pnext = PBase -> V;
  223.     else {
  224.         VBase -> Pnext = AllocVertex(0, 0, NULL, NULL);
  225.         VBase = VBase -> Pnext;
  226.         PT_COPY(VBase -> Normal, Dir);
  227.         PT_COPY(VBase -> Pt, CirclePt);
  228.     }
  229.  
  230.     PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
  231.     PT_COPY(LastCircleNrml, CircleNrml);
  232.     }
  233.  
  234.     UpdatePolyPlane(PBase, ApexPt);   /* Update base polygon plane equation. */
  235.     SET_CONVEX_POLY(PBase);               /* Mark it as convex polygon. */
  236.     PBase -> Pnext = PCone -> U.Pl.P;/* And stick it into the cone polygons. */
  237.     PCone -> U.Pl.P = PBase;
  238.  
  239.     SetObjectColor(PCone, GlblPrimColor);       /* Set its default color. */
  240.  
  241.     return PCone;
  242. }
  243.  
  244. /*****************************************************************************
  245. *   Routine to create a CONE2 geometric object defined by Pt - main base     *
  246. * 3d center point, Dir - the cone direction and length, and bases R1 and R2. *
  247. *****************************************************************************/
  248. ObjectStruct *GenCONE2Object(VectorType Pt, VectorType Dir, RealType *R1,
  249.                                 RealType *R2)
  250. {
  251.     int i, Resolution;
  252.     RealType Angle, AngleStep;
  253.     PointType LastCirclePt, CirclePt, ApexPt, LastApexPt1, ApexPt1;
  254.     NormalType LastCircleNrml, CircleNrml;
  255.     VectorType InvDir;
  256.     MatrixType Mat1, Mat2;
  257.     VertexStruct *VBase1, *VBase2, *PVertex;
  258.     PolygonStruct *PBase1, *PBase2;
  259.     ObjectStruct *PCone;
  260.  
  261.     Resolution = GetResolution(TRUE);     /* Get refinement factor of object. */
  262.  
  263.     PT_COPY(ApexPt, Pt);           /* Find the apex point: Pt + Dir. */
  264.     PT_ADD(ApexPt, ApexPt, Dir);
  265.     PT_NORMALIZE(Dir);
  266.     PT_COPY(InvDir, Dir);
  267.     PT_SCALE(InvDir, -1.0);
  268.  
  269.     GenTransformMatrix(Mat1, Pt, Dir, *R1);   /* Transform from unit circle. */
  270.     GenTransformMatrix(Mat2, ApexPt, Dir, *R2);
  271.  
  272.     PCone = GenPolyObject("", NULL, NULL);   /* Gen. the CONE object itself: */
  273.     /* Also allocate the base polygon header with first vertex on it: */
  274.     PBase1 = AllocPolygon(0, 0, VBase1 = AllocVertex(0, 0, NULL, NULL), NULL);
  275.     PBase2 = AllocPolygon(0, 0, VBase2 = AllocVertex(0, 0, NULL, NULL), NULL);
  276.  
  277.     /* First point is allways at Angle = 0. */
  278.     LastCirclePt[0] = LastApexPt1[0] = 1.0;
  279.     LastCirclePt[1] = LastApexPt1[1] = 0.0;
  280.     LastCirclePt[2] = LastApexPt1[2] = 0.0;
  281.     MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat1);
  282.     MatMultVecby4by4(LastApexPt1, LastApexPt1, Mat2);
  283.  
  284.     UpdateVertexNormal(LastCircleNrml, LastCirclePt, Pt, TRUE, ApexPt);
  285.  
  286.     PT_COPY(VBase1 -> Pt, LastCirclePt);/* Update first pt in base1 polygon. */
  287.     PT_COPY(VBase1 -> Normal, Dir);
  288.     PT_COPY(VBase2 -> Pt, LastApexPt1); /* Update first pt in base2 polygon. */
  289.     PT_COPY(VBase2 -> Normal, InvDir);
  290.  
  291.     AngleStep = M_PI * 2 / Resolution;
  292.  
  293.     for (i = 1; i <= Resolution; i++) {          /* Pass the whole base circle. */
  294.     Angle = AngleStep * i;             /* Prevent from additive error. */
  295.  
  296.     CirclePt[0] = ApexPt1[0] = cos(Angle);
  297.     CirclePt[1] = ApexPt1[1] = sin(Angle);
  298.     CirclePt[2] = ApexPt1[2] = 0.0;
  299.     MatMultVecby4by4(CirclePt, CirclePt, Mat1);
  300.     MatMultVecby4by4(ApexPt1, ApexPt1, Mat2);
  301.  
  302.      UpdateVertexNormal(CircleNrml, CirclePt, Pt, TRUE, ApexPt);
  303.  
  304.     PCone -> U.Pl.P = GenPolygon4Vrtx(LastCirclePt, LastApexPt1, ApexPt1,
  305.                       CirclePt, Pt, PCone -> U.Pl.P);
  306.  
  307.     /* Update the normals for this cone side polygon vertices: */
  308.     PVertex = PCone -> U.Pl.P -> V;
  309.     PT_COPY(PVertex -> Normal, LastCircleNrml);
  310.     PVertex = PVertex -> Pnext;
  311.     PT_COPY(PVertex -> Normal, LastCircleNrml );
  312.     PVertex = PVertex -> Pnext;
  313.     PT_COPY(PVertex -> Normal, CircleNrml);
  314.     PVertex = PVertex -> Pnext;
  315.     PT_COPY(PVertex -> Normal, CircleNrml);
  316.  
  317.     /* And add these vertices to base polygons: */
  318.     if (i == Resolution) {           /* Its last point - make it circular. */
  319.         VBase1 -> Pnext = PBase1 -> V;
  320.         VBase2 -> Pnext = PBase2 -> V;
  321.     }
  322.     else {
  323.         VBase1 -> Pnext = AllocVertex(0, 0, NULL, NULL);
  324.         VBase1 = VBase1 -> Pnext;
  325.         PT_COPY(VBase1 -> Pt, CirclePt);
  326.         PT_COPY(VBase1 -> Normal, Dir);
  327.         VBase2 -> Pnext = AllocVertex(0, 0, NULL, NULL);
  328.         VBase2 = VBase2 -> Pnext;
  329.         PT_COPY(VBase2 -> Pt, ApexPt1);
  330.         PT_COPY(VBase2 -> Normal, InvDir);
  331.     }
  332.  
  333.     PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
  334.     PT_COPY(LastApexPt1, ApexPt1);
  335.     PT_COPY(LastCircleNrml, CircleNrml);
  336.     }
  337.  
  338.     UpdatePolyPlane(PBase1, ApexPt);  /* Update base polygon plane equation. */
  339.     SET_CONVEX_POLY(PBase1);               /* Mark it as convex polygon. */
  340.     UpdatePolyPlane(PBase2, Pt);      /* Update base polygon plane equation. */
  341.     SET_CONVEX_POLY(PBase2);               /* Mark it as convex polygon. */
  342.  
  343.     PBase1 -> Pnext = PCone -> U.Pl.P;  /* And stick into the cone polygons. */
  344.     PCone -> U.Pl.P = PBase1;
  345.     PBase2 -> Pnext = PCone -> U.Pl.P;
  346.     PCone -> U.Pl.P = PBase2;
  347.  
  348.     SetObjectColor(PCone, GlblPrimColor);       /* Set its default color. */
  349.  
  350.     return PCone;
  351. }
  352.  
  353. /*****************************************************************************
  354. *   Routine to create a CYLINder geometric object defined by Pt - the base   *
  355. * 3d center point, Dir - the cylin direction and length, and base radius R.  *
  356. *   The second base is defined from first one by translating it by vector    *
  357. * Dir, and its points are prefixed with T.                     *
  358. *****************************************************************************/
  359. ObjectStruct *GenCYLINObject(VectorType Pt, VectorType Dir, RealType *R)
  360. {
  361.     int i, Resolution;
  362.     RealType Angle, AngleStep;
  363.     PointType LastCirclePt, CirclePt, TLastCirclePt, TCirclePt, TPt, Dummy;
  364.     VectorType ForwardDir, BackwardDir;
  365.     NormalType LastCircleNrml, CircleNrml;
  366.     MatrixType Mat;
  367.     VertexStruct *VBase1, *VBase2, *PVertex;
  368.     PolygonStruct *PBase1, *PBase2;
  369.     ObjectStruct *PCylin;
  370.  
  371.     Resolution = GetResolution(TRUE);     /* Get refinement factor of object. */
  372.  
  373.     GenTransformMatrix(Mat, Pt, Dir, *R);     /* Transform from unit circle. */
  374.  
  375.     PCylin = GenPolyObject("", NULL, NULL); /* Gen. the CYLIN object itself: */
  376.     /* Also allocate the bases polygon header with first vertex on it: */
  377.     PBase1 = AllocPolygon(0, 0, VBase1 = AllocVertex(0, 0, NULL, NULL), NULL);
  378.     PBase2 = AllocPolygon(0, 0, VBase2 = AllocVertex(0, 0, NULL, NULL), NULL);
  379.  
  380.     PT_ADD(TPt, Pt, Dir);           /* Translated circle center (by Dir). */
  381.  
  382.     /* Prepare the normal directions for the two bases: */
  383.     PT_COPY(ForwardDir, Dir);
  384.     PT_NORMALIZE(ForwardDir);
  385.     PT_COPY(BackwardDir, ForwardDir);
  386.     PT_SCALE(BackwardDir, -1.0);
  387.  
  388.     LastCirclePt[0] = 1.0;        /* First point is allways Angle = 0. */
  389.     LastCirclePt[1] = 0.0;
  390.     LastCirclePt[2] = 0.0;
  391.     MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);
  392.  
  393.     UpdateVertexNormal(LastCircleNrml, LastCirclePt, Pt, FALSE, Dummy);
  394.  
  395.     PT_COPY(VBase1 -> Pt, LastCirclePt);/* Update first pt in base1 polygon. */
  396.     PT_COPY(VBase1 -> Normal, ForwardDir);
  397.     PT_ADD(TLastCirclePt, LastCirclePt, Dir); /* Translated circle (by Dir). */
  398.     PT_COPY(VBase2 -> Pt, TLastCirclePt);/* Update first pt in base2 polygon.*/
  399.     PT_COPY(VBase2 -> Normal, BackwardDir);
  400.  
  401.     AngleStep = M_PI * 2 / Resolution;
  402.  
  403.     for (i = 1; i <= Resolution; i++) {          /* Pass the whole base circle. */
  404.     Angle = AngleStep * i;             /* Prevent from additive error. */
  405.  
  406.     CirclePt[0] = cos(Angle);
  407.     CirclePt[1] = sin(Angle);
  408.     CirclePt[2] = 0.0;
  409.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  410.  
  411.     UpdateVertexNormal(CircleNrml, CirclePt, Pt, FALSE, Dummy);
  412.  
  413.     PT_ADD(TCirclePt, CirclePt, Dir);     /* Translated circle (by Dir). */
  414.  
  415.     PCylin -> U.Pl.P = GenPolygon4Vrtx(TLastCirclePt, TCirclePt, CirclePt,
  416.                      LastCirclePt, Pt, PCylin -> U.Pl.P);
  417.     /* Update the normals for this cylinder side polygon vertices: */
  418.     PVertex = PCylin -> U.Pl.P -> V;
  419.     PT_COPY(PVertex -> Normal, LastCircleNrml);
  420.     PVertex = PVertex -> Pnext;
  421.     PT_COPY(PVertex -> Normal, CircleNrml);
  422.     PVertex = PVertex -> Pnext;
  423.     PT_COPY(PVertex -> Normal, CircleNrml);
  424.     PVertex = PVertex -> Pnext;
  425.     PT_COPY(PVertex -> Normal, LastCircleNrml);
  426.  
  427.     /* And add this vertices to the two cylinder bases:              */
  428.     /* Note Base1 is build forward, while Base2 is build backward so it  */
  429.     /* will be consistent - cross product of 2 consecutive edges will    */
  430.     /* point into the model.                         */
  431.     if (i == Resolution) {           /* Its last point - make it circular. */
  432.         VBase1 -> Pnext = PBase1 -> V;
  433.         VBase2 -> Pnext = PBase2 -> V;
  434.     }
  435.     else {
  436.         VBase1 -> Pnext = AllocVertex(0, 0, NULL, NULL);
  437.         VBase1 = VBase1 -> Pnext;
  438.         PT_COPY(VBase1 -> Pt, CirclePt);
  439.         PT_COPY(VBase1 -> Normal, ForwardDir);
  440.         PBase2 -> V = AllocVertex(0, 0, NULL, PBase2 -> V);
  441.         PT_COPY(PBase2 -> V -> Pt, TCirclePt);
  442.         PT_COPY(PBase2 -> V -> Normal, BackwardDir);
  443.     }
  444.  
  445.     PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
  446.     PT_COPY(TLastCirclePt, TCirclePt);
  447.     PT_COPY(LastCircleNrml, CircleNrml);
  448.     }
  449.  
  450.     UpdatePolyPlane(PBase1, TPt);     /* Update base polygon plane equation. */
  451.     SET_CONVEX_POLY(PBase1);               /* Mark it as convex polygon. */
  452.     PBase1 -> Pnext = PCylin -> U.Pl.P; /* And stick it into cylin polygons. */
  453.     PCylin -> U.Pl.P = PBase1;
  454.     UpdatePolyPlane(PBase2, Pt);      /* Update base polygon plane equation. */
  455.     SET_CONVEX_POLY(PBase2);               /* Mark it as convex polygon. */
  456.     PBase2 -> Pnext = PCylin -> U.Pl.P; /* And stick it into cylin polygons. */
  457.     PCylin -> U.Pl.P = PBase2;
  458.  
  459.     SetObjectColor(PCylin, GlblPrimColor);       /* Set its default color. */
  460.  
  461.     return PCylin;
  462. }
  463.  
  464. /*****************************************************************************
  465. *   Routine to create a SPHERE geometric object defined by Center - sphere   *
  466. * 3d center point, and R its radius.                         *
  467. *   Note polygons on the poles are triangles, while others are rectangles    *
  468. *   Teta is horizontal circle angle, Fee is the vertical (spherical coords.) *
  469. *   The vertical axes here is assumed to be the Z axes.                 *
  470. *****************************************************************************/
  471. ObjectStruct *GenSPHEREObject(VectorType Center, RealType *R)
  472. {
  473.     int i, j, k, Resolution;
  474.     RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
  475.     CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
  476.     PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt, Dummy;
  477.     VertexStruct *PVertex;
  478.     ObjectStruct *PSphere;
  479.  
  480.     Resolution = GetResolution(TRUE);     /* Get refinement factor of object. */
  481.  
  482.     PSphere = GenPolyObject("", NULL, NULL);/* Gen the SPHERE object itself: */
  483.  
  484.     TetaAngleStep = M_PI * 2.0 / Resolution;         /* Runs from 0 to 2*PI. */
  485.     FeeAngleStep = M_PI * 2.0 / Resolution;    /* Runs from -PI/2 yo +PI/2. */
  486.  
  487.     /* Generate the lowest (south pole) triangular polygons: */
  488.     FeeAngle = (-M_PI/2.0) + FeeAngleStep; /* First circle above south pole. */
  489.     CosFeeAngle1 = cos(FeeAngle) * (*R);
  490.     SinFeeAngle1 = sin(FeeAngle) * (*R);
  491.     PT_COPY(LastCirclePt, Center);        /* Calculate the south pole. */
  492.     LastCirclePt[2] -= (*R);
  493.     PT_COPY(CircleLastPt, Center);    /* Calc. last point on current circle. */
  494.     CircleLastPt[0] += CosFeeAngle1;
  495.     CircleLastPt[2] += SinFeeAngle1;
  496.  
  497.     for (i = 1; i <= Resolution; i++) {   /* Pass whole (horizontal) circle. */
  498.     TetaAngle = TetaAngleStep * i;         /* Prevent from additive error. */
  499.  
  500.     PT_COPY(CirclePt, Center); /* Calc. current point on current circle. */
  501.     CirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
  502.     CirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
  503.     CirclePt[2] += SinFeeAngle1;
  504.  
  505.     PSphere -> U.Pl.P = GenPolygon3Vrtx(LastCirclePt, CircleLastPt,
  506.                     CirclePt, Center, PSphere -> U.Pl.P);
  507.     /* Update normals: */
  508.     for (j = 0, PVertex = PSphere -> U.Pl.P -> V;
  509.          j < 3;
  510.          j++, PVertex = PVertex -> Pnext)
  511.         UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, Center,
  512.                                 FALSE, Dummy);
  513.  
  514.     PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
  515.     }
  516.  
  517.     /* Generate the middle rectangular polygons: */
  518.     for (i = 1; i < Resolution/2-1; i++) { /* For all horizontal circles do. */
  519.     FeeAngle = (-M_PI/2.0) + FeeAngleStep * i;
  520.     CosFeeAngle1 = cos(FeeAngle) * (*R);
  521.     SinFeeAngle1 = sin(FeeAngle) * (*R);
  522.     FeeAngle = (-M_PI/2.0) + FeeAngleStep * (i + 1);
  523.     CosFeeAngle2 = cos(FeeAngle) * (*R);
  524.     SinFeeAngle2 = sin(FeeAngle) * (*R);
  525.     PT_COPY(CircleLastPt, Center);/* Calc. last point on current circle. */
  526.     CircleLastPt[0] += CosFeeAngle2;
  527.     CircleLastPt[2] += SinFeeAngle2;
  528.     PT_COPY(LastCircleLastPt, Center);/* Calc. last point on last circle.*/
  529.     LastCircleLastPt[0] += CosFeeAngle1;
  530.     LastCircleLastPt[2] += SinFeeAngle1;
  531.  
  532.     for (j = 1; j <= Resolution; j++) {/* Pass whole (horizontal) circle.*/
  533.         TetaAngle = TetaAngleStep * j;   /* Prevent from additive error. */
  534.  
  535.         PT_COPY(CirclePt, Center);/* Calc. current pt on current circle. */
  536.         CirclePt[0] += cos(TetaAngle) * CosFeeAngle2;
  537.         CirclePt[1] += sin(TetaAngle) * CosFeeAngle2;
  538.         CirclePt[2] += SinFeeAngle2;
  539.         PT_COPY(LastCirclePt, Center);/* Calc. current pt on last circle.*/
  540.         LastCirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
  541.         LastCirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
  542.         LastCirclePt[2] += SinFeeAngle1;
  543.  
  544.         PSphere -> U.Pl.P = GenPolygon4Vrtx(LastCirclePt, LastCircleLastPt,
  545.             CircleLastPt, CirclePt, Center, PSphere -> U.Pl.P);
  546.         /* Update normals: */
  547.         for (k = 0, PVertex = PSphere -> U.Pl.P -> V;
  548.          k < 4;
  549.          k++, PVertex = PVertex -> Pnext)
  550.         UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, Center,
  551.                                 FALSE, Dummy);
  552.  
  553.         PT_COPY(CircleLastPt, CirclePt);          /* Save pt in last pt. */
  554.         PT_COPY(LastCircleLastPt, LastCirclePt);
  555.     }
  556.     }
  557.  
  558.     /* Generate the upper most (north pole) triangular polygons: */
  559.     FeeAngle = (M_PI/2.0) - FeeAngleStep;  /* First circle below north pole. */
  560.     CosFeeAngle1 = cos(FeeAngle) * (*R);
  561.     SinFeeAngle1 = sin(FeeAngle) * (*R);
  562.     PT_COPY(LastCirclePt, Center);        /* Calculate the north pole. */
  563.     LastCirclePt[2] += (*R);
  564.     PT_COPY(CircleLastPt, Center);    /* Calc. last point on current circle. */
  565.     CircleLastPt[0] += CosFeeAngle1;
  566.     CircleLastPt[2] += SinFeeAngle1;
  567.  
  568.     for (i = 1; i <= Resolution; i++) {   /* Pass whole (horizontal) circle. */
  569.     TetaAngle = TetaAngleStep * i;         /* Prevent from additive error. */
  570.  
  571.     PT_COPY(CirclePt, Center); /* Calc. current point on current circle. */
  572.     CirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
  573.     CirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
  574.     CirclePt[2] += SinFeeAngle1;
  575.  
  576.     PSphere -> U.Pl.P = GenPolygon3Vrtx(LastCirclePt, CirclePt, CircleLastPt,
  577.                         Center, PSphere -> U.Pl.P);
  578.  
  579.     /* Update normals: */
  580.     for (j = 0, PVertex = PSphere -> U.Pl.P -> V;
  581.          j < 3;
  582.          j++, PVertex = PVertex -> Pnext)
  583.         UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, Center,
  584.                                 FALSE, Dummy);
  585.  
  586.     PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
  587.     }
  588.  
  589.     SetObjectColor(PSphere, GlblPrimColor);       /* Set its default color. */
  590.  
  591.     return PSphere;
  592. }
  593.  
  594. /*****************************************************************************
  595. *   Routine to create a TORUS geometric object defined by Center - torus 3d  *
  596. * center point, the main torus plane normal Normal, major radius Rmajor and  *
  597. * minor radius Rminor (Tube radius).                         *
  598. *   Teta runs on the major circle, Fee on the minor one (Before Mat trans.): *
  599. * X = (Rmajor + Rminor * cos(Fee)) * cos(Teta)                     *
  600. * Y = (Rmajor + Rminor * cos(Fee)) * sin(Teta)                     *
  601. * Z = Rminor * sin(Fee)                                 *
  602. *****************************************************************************/
  603. ObjectStruct *GenTORUSObject(VectorType Center, VectorType Normal,
  604.                     RealType *Rmajor, RealType *Rminor)
  605. {
  606.     int i, j, Resolution;
  607.     RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
  608.     CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
  609.     PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt,
  610.     LastInPt, InPt, Dummy;
  611.     MatrixType Mat;
  612.     VertexStruct *PVertex;
  613.     ObjectStruct *PTorus;
  614.  
  615.     Resolution = GetResolution(TRUE);     /* Get refinement factor of object. */
  616.  
  617.     GenTransformMatrix(Mat, Center, Normal, 1.0); /* Trans from unit circle. */
  618.  
  619.     PTorus = GenPolyObject("", NULL, NULL); /* Gen. the Torus object itself: */
  620.  
  621.     TetaAngleStep = M_PI * 2.0 / Resolution;         /* Runs from 0 to 2*PI. */
  622.     FeeAngleStep = M_PI * 2.0 / Resolution;         /* Runs from 0 to 2*PI. */
  623.  
  624.     for (i = 1; i <= Resolution; i++) {
  625.     FeeAngle = FeeAngleStep * (i - 1);
  626.     CosFeeAngle1 = cos(FeeAngle) * (*Rminor);
  627.     SinFeeAngle1 = sin(FeeAngle) * (*Rminor);
  628.     FeeAngle = FeeAngleStep * i;
  629.     CosFeeAngle2 = cos(FeeAngle) * (*Rminor);
  630.     SinFeeAngle2 = sin(FeeAngle) * (*Rminor);
  631.     LastCircleLastPt[0] = (*Rmajor) + CosFeeAngle1;
  632.     LastCircleLastPt[1] = 0.0;
  633.     LastCircleLastPt[2] = SinFeeAngle1;
  634.     MatMultVecby4by4(LastCircleLastPt, LastCircleLastPt, Mat);
  635.     LastCirclePt[0] = (*Rmajor) + CosFeeAngle2;
  636.     LastCirclePt[1] = 0.0;
  637.     LastCirclePt[2] = SinFeeAngle2;
  638.     MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);
  639.     /* Point inside the object relative to this polygon: */
  640.     LastInPt[0] = (*Rmajor);
  641.     LastInPt[1] = 0.0;
  642.     LastInPt[2] = 0.0;
  643.     MatMultVecby4by4(LastInPt, LastInPt, Mat);
  644.  
  645.     for (j = 1; j <= Resolution; j++) {
  646.         TetaAngle = TetaAngleStep * j;   /* Prevent from additive error. */
  647.  
  648.         CircleLastPt[0] = ((*Rmajor) + CosFeeAngle1) * cos(TetaAngle);
  649.         CircleLastPt[1] = ((*Rmajor) + CosFeeAngle1) * sin(TetaAngle);
  650.         CircleLastPt[2] = SinFeeAngle1;
  651.         MatMultVecby4by4(CircleLastPt, CircleLastPt, Mat);
  652.         CirclePt[0] = ((*Rmajor) + CosFeeAngle2) * cos(TetaAngle);
  653.         CirclePt[1] = ((*Rmajor) + CosFeeAngle2) * sin(TetaAngle);
  654.         CirclePt[2] = SinFeeAngle2;
  655.         MatMultVecby4by4(CirclePt, CirclePt, Mat);
  656.         /* Point inside the object relative to this polygon: */
  657.         InPt[0] = (*Rmajor) * cos(TetaAngle);
  658.         InPt[1] = (*Rmajor) * sin(TetaAngle);
  659.         InPt[2] = 0.0;
  660.         MatMultVecby4by4(InPt, InPt, Mat);
  661.  
  662.         PTorus -> U.Pl.P = GenPolygon4Vrtx(CirclePt, CircleLastPt,
  663.         LastCircleLastPt, LastCirclePt, InPt, PTorus -> U.Pl.P);
  664.  
  665.         /* Update normals: */
  666.         PVertex = PTorus -> U.Pl.P -> V;
  667.         UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, InPt,
  668.                                 FALSE, Dummy);
  669.         PVertex = PVertex -> Pnext;
  670.         UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, InPt,
  671.                                 FALSE, Dummy);
  672.         PVertex = PVertex -> Pnext;
  673.         UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, LastInPt,
  674.                                 FALSE, Dummy);
  675.         PVertex = PVertex -> Pnext;
  676.         UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, LastInPt,
  677.                                 FALSE, Dummy);
  678.  
  679.         PT_COPY(LastCirclePt, CirclePt);          /* Save pt in last pt. */
  680.         PT_COPY(LastCircleLastPt, CircleLastPt);
  681.         PT_COPY(LastInPt, InPt);
  682.     }
  683.     }
  684.  
  685.     SetObjectColor(PTorus, GlblPrimColor);       /* Set its default color. */
  686.  
  687.     return PTorus;
  688. }
  689.  
  690. /*****************************************************************************
  691. *   Routine to create a PLANE geometric object defined by the normal N and   *
  692. * the translation vector T. The object is a FINITE plane (a circle of         *
  693. * RESOLUTION point in it...) and its radius is equal to R.             *
  694. *   The normal direction is assumed to point to the inside of the object.    *
  695. *****************************************************************************/
  696. ObjectStruct *GenPLANEObject(VectorType N, VectorType T, RealType *R)
  697. {
  698.     int i, Resolution;
  699.     RealType Angle, AngleStep;
  700.     PointType CirclePt;
  701.     MatrixType Mat;
  702.     VertexStruct *V;
  703.     PolygonStruct *PCirc;
  704.     ObjectStruct *PPlane;
  705.  
  706.     Resolution = GetResolution(TRUE);     /* Get refinement factor of object. */
  707.  
  708.     GenTransformMatrix(Mat, T, N, *R);          /* Transform from unit circle. */
  709.     PT_NORMALIZE(N);
  710.  
  711.     PPlane = GenPolyObject("", NULL, NULL); /* Gen. the PLANE object itself: */
  712.     PCirc = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), NULL);
  713.     PPlane -> U.Pl.P = PCirc;
  714.  
  715.     CirclePt[0] = 1.0;            /* First point is allways Angle = 0. */
  716.     CirclePt[1] = 0.0;
  717.     CirclePt[2] = 0.0;
  718.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  719.     PT_COPY(V -> Pt, CirclePt);           /* Update first pt in circle polygon. */
  720.     PT_COPY(V -> Normal, N);
  721.  
  722.     AngleStep = M_PI * 2 / Resolution;
  723.  
  724.     for (i = 1; i <= Resolution; i++) {          /* Pass the whole base circle. */
  725.     Angle = AngleStep * i;             /* Prevent from additive error. */
  726.  
  727.     CirclePt[0] = cos(Angle);
  728.     CirclePt[1] = sin(Angle);
  729.     CirclePt[2] = 0.0;
  730.  
  731.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  732.  
  733.     /* And add this vertices to the two cylinder bases: */
  734.     if (i == Resolution) {           /* Its last point - make it circular. */
  735.         V -> Pnext = PCirc -> V;
  736.     }
  737.     else {
  738.         V -> Pnext = AllocVertex(0, 0, NULL, NULL);
  739.         V = V -> Pnext;
  740.         PT_COPY(V -> Pt, CirclePt);
  741.         PT_COPY(V -> Normal, N);
  742.     }
  743.     }
  744.  
  745.     PT_ADD(CirclePt, CirclePt, N);   /* Make a point "IN" the circle object. */
  746.     UpdatePolyPlane(PCirc, CirclePt); /* Update base polygon plane equation. */
  747.     SET_CONVEX_POLY(PCirc);               /* Mark it as convex polygon. */
  748.  
  749.     SetObjectColor(PPlane, GlblPrimColor);       /* Set its default color. */
  750.  
  751.     return PPlane;
  752. }
  753.  
  754. /*****************************************************************************
  755. *   Routine to create a POLYGON directly from its specified vertices.         *
  756. *   The validity of list elements is tested to make sure they are vectors as *
  757. * nobody else checked it till now (lists can hold any object type!).         *
  758. *   No test is made to make sure all vertices are on one plane, and that no  *
  759. * two vertices are similar.                             *
  760. *****************************************************************************/
  761. ObjectStruct *GenPOLYGONObject(ObjectStruct *PObjList)
  762. {
  763.     int i, NumVertices = 0;
  764.     VertexStruct *V, *VHead, *VTail;
  765.     PolygonStruct *PPoly;
  766.     ObjectStruct *PObj, *PObjPoly;
  767.  
  768.     if (!IS_OLST_OBJ(PObjList))
  769.     FatalError("GenPOLYObject: Not object list object!");
  770.  
  771.     i = 0;
  772.     while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
  773.     if (!IS_VEC_OBJ(PObj)) {
  774.         WndwInputWindowPutStr("POLY: None vector object found in list, empty object result.");
  775.         return NULL;
  776.     }
  777.     NumVertices++;
  778.     }
  779.  
  780.     if (NumVertices < 3) {
  781.     WndwInputWindowPutStr("POLY: Less than 3 vertices, empty object result.");
  782.     return NULL;
  783.     }
  784.  
  785.     PPoly = AllocPolygon(0, 0, VHead = NULL, NULL);
  786.     i = 0;
  787.     while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
  788.     V = AllocVertex(0, 0, NULL, NULL);
  789.     PT_COPY(V -> Pt, PObj -> U.Vec);
  790.     if (VHead == NULL) {
  791.         PPoly -> V = VHead = VTail = V;
  792.     }
  793.     else {
  794.         VTail -> Pnext = V;
  795.         VTail = V;
  796.     }
  797.     }
  798.     VTail -> Pnext = VHead;              /* Close the vertex list loop. */
  799.  
  800.     /* Update polygon plane equation and vertices normals. */
  801.     CGPlaneFrom3Points(PPoly -> Plane, VHead -> Pt, VHead -> Pnext -> Pt,
  802.                         VHead -> Pnext -> Pnext -> Pt);
  803.     V = VHead;
  804.     do {
  805.     PT_COPY(V -> Normal, PPoly -> Plane);
  806.     V = V -> Pnext;
  807.     }
  808.     while (V != VHead);
  809.  
  810.     PObjPoly = GenPolyObject("", NULL, NULL);/* Gen. the POLY object itself: */
  811.     PObjPoly -> U.Pl.P = PPoly;
  812.     SetObjectColor(PObjPoly, GlblPrimColor);       /* Set its default color. */
  813.  
  814.     return PObjPoly;
  815. }
  816.  
  817. /*****************************************************************************
  818. *   Routine to create an OBJECT directly from set of sspecified polygons.    *
  819. *   No test is made for the validity of the model in any sense.             *
  820. *****************************************************************************/
  821. ObjectStruct *GenObjectFromPolyList(ObjectStruct *PObjList)
  822. {
  823.     int i;
  824.     PolygonStruct *PPoly,
  825.     *PTail = NULL;
  826.     ObjectStruct *PObj, *PObjPoly;
  827.  
  828.     if (!IS_OLST_OBJ(PObjList))
  829.     FatalError("GenObjectFromPolyList: Not object list object!");
  830.  
  831.     i = 0;
  832.     while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
  833.     if (!IS_POLY_OBJ(PObj)) {
  834.         WndwInputWindowPutStr("POLYTOOBJ: None polygon object found in list, empty object result.");
  835.         return NULL;
  836.     }
  837.     }
  838.  
  839.     PObjPoly = GenPolyObject("", NULL, NULL);/* Gen. the POLY object itself: */
  840.     SetObjectColor(PObjPoly, GlblPrimColor);       /* Set its default color. */
  841.     i = 0;
  842.     while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
  843.     PPoly = CopyPolygonList(PObj -> U.Pl.P);
  844.  
  845.     if (PTail == NULL) {
  846.         PObjPoly -> U.Pl.P = PPoly;
  847.     }
  848.     else {
  849.         PTail -> Pnext = PPoly;
  850.     }
  851.     for (PTail = PPoly; PTail -> Pnext != NULL; PTail = PTail -> Pnext);
  852.     }
  853.  
  854.     return PObjPoly;
  855. }
  856.  
  857. /*****************************************************************************
  858. *   Routine to a cross section polygon, by interactively allowing the user   *
  859. * to create it. This polygon is returned as a one polygon object. This       *
  860. * polygon is mainly for Surface of Revolution (SURFREV) and extrusion        *
  861. * (EXTRUDE) operations, although it may be used as an open object by itself. *
  862. *****************************************************************************/
  863. ObjectStruct *GenCROSSECObject(ObjectStruct *PObj)
  864. {
  865.     if (PObj && !IS_POLY_OBJ(PObj))
  866.     FatalError("CrossSec: operation on non polygonal object");
  867.  
  868.     WndwInputWindowPutStr("GenCrossSecObject not implemented");
  869.  
  870.     return NULL;
  871. }
  872.  
  873. /*****************************************************************************
  874. *   Routine to a create surface of revolution by rotating the given cross    *
  875. * section along Z axes.                                 *
  876. * If input is a polyline/gon, it must not be perpendicular to Z.         *
  877. *****************************************************************************/
  878. ObjectStruct *GenSURFREVObject(ObjectStruct *Cross)
  879. {
  880.     int Resolution, i, j;
  881.     RealType XYSize;
  882.     MatrixType Mat;               /* Rotation Matrix around Z axes. */
  883.     VertexStruct *V, *V1, *V1Head, *V2, *V2Head, *VIn, *VInHead;
  884.     PolygonStruct *Pl1, *Pl2, *PlIn, *PlNew = NULL;
  885.     ObjectStruct *PSurfRev;
  886.  
  887.     if (!IS_POLY_OBJ(Cross) && !IS_CRV_OBJ(Cross)) {
  888.     WndwInputWindowPutStr("SurfRev: cross section is not poly/crv. Empty object result");
  889.     return NULL;
  890.     }
  891.  
  892.     if (IS_POLY_OBJ(Cross)) {
  893.     if (APX_EQ(Cross -> U.Pl.P -> Plane[0], 0.0) &&
  894.         APX_EQ(Cross -> U.Pl.P -> Plane[1], 0.0)) {
  895.         WndwInputWindowPutStr("SurfRev: cross-section perpendicular to Z. Empty object result");
  896.         return NULL;
  897.     }
  898.  
  899.     Pl1 = AllocPolygon(0, 0, V1Head = CopyVList(Cross -> U.Pl.P -> V), NULL);
  900.     PLANE_COPY(Pl1 -> Plane, Cross -> U.Pl.P -> Plane);
  901.     Pl2 = AllocPolygon(0, 0, V2Head = CopyVList(Cross -> U.Pl.P -> V), NULL);
  902.     PLANE_COPY(Pl2 -> Plane, Cross -> U.Pl.P -> Plane);
  903.     PlIn = GenInsidePoly(Pl1);
  904.     VInHead = PlIn -> V;
  905.     Resolution = GetResolution(TRUE);/* Get refinement factor of object. */
  906.     MatGenMatRotZ1(2.0 * M_PI / Resolution, Mat);
  907.  
  908.     for (i = 0; i < Resolution; i++)
  909.     {
  910.         V2 = V2Head;
  911.         do {
  912.         MatMultVecby4by4(V2 -> Pt, V2 -> Pt , Mat);
  913.         V2 = V2 -> Pnext;
  914.         }
  915.         while (V2 != NULL && V2 != V2Head);
  916.  
  917.         V1 = V1Head;
  918.         if (i < Resolution - 1) /* If this is last loop use the original */
  919.             V2 = V2Head; /* poly as we might accumulate error during the */
  920.         else            /* transformations along the circle. */
  921.         V2 = Cross -> U.Pl.P -> V;
  922.         VIn = VInHead;
  923.  
  924.         do {
  925.         PlNew = GenPolygon4Vrtx(V1 -> Pt, V1 -> Pnext -> Pt,
  926.                         V2 -> Pnext -> Pt, V2 -> Pt,
  927.                         VIn -> Pt, PlNew);
  928.  
  929.             /* Update normals: */
  930.             for (j = 0, V = PlNew -> V; j < 4; j++, V = V -> Pnext) {
  931.             V -> Normal[0] = V -> Pt[0];
  932.             V -> Normal[1] = V -> Pt[1];
  933.             V -> Normal[2] = 0.0;
  934.             /* Make sure normal does not point in opposite direction.*/
  935.             if (DOT_PROD(V -> Normal, PlNew -> Plane) < 0.0)
  936.                 PT_SCALE(V -> Normal, -1.0);
  937.  
  938.             /* Since Z normal component should be fixed for all normals: */
  939.             V -> Normal[2] = PlNew -> Plane[2];
  940.             XYSize = APX_EQ(ABS(PlNew -> Plane[2]), 1.0) ?
  941.                     0.0 : 1 - SQR(PlNew -> Plane[2]);
  942.             XYSize = sqrt(XYSize / (SQR(V -> Pt[0]) + SQR(V -> Pt[1])));
  943.             V -> Normal[0] *= XYSize;
  944.             V -> Normal[1] *= XYSize;
  945.             }
  946.  
  947.             VIn = VIn -> Pnext;
  948.             V1 = V1 -> Pnext;
  949.             V2 = V2 -> Pnext;
  950.         }
  951.         while (V1 -> Pnext != NULL && V1 != V1Head);
  952.  
  953.         V1 = V1Head;
  954.         do {
  955.             MatMultVecby4by4(V1 -> Pt, V1 -> Pt , Mat);
  956.             V1 = V1 -> Pnext;
  957.         }
  958.         while (V1 != NULL && V1 != V1Head);
  959.         VIn = VInHead;
  960.         do {
  961.             MatMultVecby4by4(VIn -> Pt, VIn -> Pt , Mat);
  962.             VIn = VIn -> Pnext;
  963.         }
  964.         while (VIn != NULL && VIn != VInHead);
  965.         }
  966.  
  967.         MyFree((char *) PlIn, ALLOC_POLYGON);
  968.         MyFree((char *) Pl1, ALLOC_POLYGON);
  969.         MyFree((char *) Pl2, ALLOC_POLYGON);
  970.  
  971.         PSurfRev = GenPolyObject("", NULL, NULL);
  972.     PSurfRev -> U.Pl.P = PlNew;
  973.     SetObjectColor(PSurfRev, GlblPrimColor);   /* Set its default color. */
  974.  
  975.         return PSurfRev;
  976.     }
  977.     else if (IS_CRV_OBJ(Cross)) {
  978.     if (CAGD_NUM_OF_PT_COORD(Cross->U.Crv.Crv -> PType) < 3) {
  979.         WndwInputWindowPutStr("SurfRev: cross-section perpendicular to Z. Empty object result");
  980.         return NULL;
  981.     }
  982.  
  983.         PSurfRev = GenSrfObject("", CagdSurfaceRev(Cross->U.Crv.Crv), NULL);
  984.     SetObjectColor(PSurfRev, GlblPrimColor);   /* Set its default color. */
  985.     return PSurfRev;
  986.     }
  987.     else
  988.         return NULL;
  989. }
  990.  
  991. /*****************************************************************************
  992. *   Routine to a create surface of extrusion out of the given cross section  *
  993. * and the given direction. A NULL object will be returned, if the direction  *
  994. * vector is in the cross section plane.                         *
  995. *****************************************************************************/
  996. ObjectStruct *GenEXTRUDEObject(ObjectStruct *Cross, VectorType Dir)
  997. {
  998.     int i;
  999.     RealType R;
  1000.     CagdVecStruct CagdDir;
  1001.     VertexStruct *V1, *V1Head, *V2, *VIn;
  1002.     PolygonStruct *PBase1, *PBase2, *Pl, *PlIn;
  1003.     ObjectStruct *PExtrude;
  1004.  
  1005.     if (!IS_POLY_OBJ(Cross) && !IS_CRV_OBJ(Cross)) {
  1006.     WndwInputWindowPutStr("Extrude: cross section is not poly/crv. Empty object result");
  1007.     return NULL;
  1008.     }
  1009.  
  1010.     if (IS_POLY_OBJ(Cross)) {
  1011.     R = DOT_PROD(Cross -> U.Pl.P -> Plane, Dir);
  1012.     if (APX_EQ(R, 0.0)) {
  1013.         WndwInputWindowPutStr("Extrusion direction in cross-section plane. Empty object result");
  1014.         return NULL;
  1015.     }
  1016.  
  1017.     /* Prepare two bases (update their plane normal to point INDISE): */
  1018.     PBase1 = AllocPolygon(0, 0, CopyVList(Cross -> U.Pl.P -> V), NULL);
  1019.     Pl = PBase2 = AllocPolygon(0, 0, CopyVList(Cross -> U.Pl.P -> V), PBase1);
  1020.     V1 = V1Head = PBase2 -> V;
  1021.     do {
  1022.         PT_ADD(V1 -> Pt, Dir, V1 -> Pt);
  1023.         V1 = V1 -> Pnext;
  1024.     }
  1025.     while (V1 != NULL && V1 != V1Head);
  1026.     if (R > 0.0) {
  1027.         PLANE_COPY(PBase1 -> Plane, Cross -> U.Pl.P -> Plane);
  1028.         for (i = 0; i < 3; i++)
  1029.         PBase2 -> Plane[i] = (-Cross -> U.Pl.P -> Plane[i]);
  1030.         PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane, PBase2 -> V -> Pt));
  1031.     }
  1032.     else {
  1033.         for (i = 0; i < 4; i++)
  1034.         PBase1 -> Plane[i] = (-Cross -> U.Pl.P -> Plane[i]);
  1035.         PLANE_COPY(PBase2 -> Plane, Cross -> U.Pl.P -> Plane);
  1036.         PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane, PBase2 -> V -> Pt));
  1037.     }
  1038.  
  1039.     /* Now generate all the 4 corner polygon between the two bases: */
  1040.     V1 = V1Head = PBase1 -> V;
  1041.     V2 = PBase2 -> V;
  1042.     PlIn = GenInsidePoly(PBase1);
  1043.     VIn = PlIn -> V;
  1044.     do {
  1045.         Pl = GenPolygon4Vrtx(V1 -> Pt, V1 -> Pnext -> Pt,
  1046.                      V2 -> Pnext -> Pt, V2 -> Pt, VIn -> Pt, Pl);
  1047.         VIn = VIn -> Pnext;
  1048.         V1 = V1 -> Pnext;
  1049.         V2 = V2 -> Pnext;
  1050.     }
  1051.     while (V1 -> Pnext != NULL && V1 != V1Head);
  1052.  
  1053.     MyFree((char *) PlIn, ALLOC_POLYGON);
  1054.  
  1055.     PExtrude = GenPolyObject("", NULL, NULL);
  1056.     PExtrude -> U.Pl.P = Pl;
  1057.     SetObjectColor(PExtrude, GlblPrimColor);   /* Set its default color. */
  1058.  
  1059.     /* Update all the polygon vertices normals. */
  1060.     for (Pl = PExtrude -> U.Pl.P; Pl != NULL; Pl = Pl -> Pnext) {
  1061.         V1 = V1Head = Pl -> V;
  1062.             do {
  1063.                 PT_COPY(V1 -> Normal, Pl -> Plane);
  1064.                 V1 = V1 -> Pnext;
  1065.             }
  1066.             while (V1 != NULL && V1 != V1Head);
  1067.     }
  1068.  
  1069.     return PExtrude;
  1070.     }
  1071.     else if (IS_CRV_OBJ(Cross)) {
  1072.         for (i = 0; i < 3; i++) CagdDir.Vec[i] = Dir[i];
  1073.         PExtrude = GenSrfObject("", CagdExtrudeSrf(Cross->U.Crv.Crv, &CagdDir),
  1074.                                          NULL);
  1075.     SetObjectColor(PExtrude, GlblPrimColor);   /* Set its default color. */
  1076.         return PExtrude;
  1077.     }
  1078.     else
  1079.         return NULL;
  1080. }
  1081.  
  1082. /*****************************************************************************
  1083. *   Routine to create a pseudo polygon out of a given polygon such that each *
  1084. * vertex Vi is in the inside side of the corresponding edge ViVi+1 in the    *
  1085. * given polygon. Used in polygon generation for EXTRUDE/SURFREV operations.  *
  1086. *****************************************************************************/
  1087. static PolygonStruct *GenInsidePoly(PolygonStruct *Pl)
  1088. {
  1089.     int Axes;
  1090.     RealType Dx, Dy;
  1091.     PointType Pt;
  1092.     MatrixType Mat;
  1093.     PolygonStruct *PlIn;
  1094.     VertexStruct *VHead, *V, *Vnext, *VInHead, *VIn = NULL;
  1095.  
  1096.     PlIn = AllocPolygon(0, 0, VInHead = NULL, NULL);
  1097.  
  1098.     /* Generate transformation matrix to bring polygon to a XY parallel      */
  1099.     /* plane, and transform a copy of the polygon to that plane.         */
  1100.     GenRotateMatrix(Mat, Pl -> Plane);
  1101.     VHead = V = CopyVList(Pl -> V);      /* We dont want to modify original! */
  1102.     Pl = AllocPolygon(0, 0, VHead, NULL);
  1103.     do {
  1104.     MatMultVecby4by4(V -> Pt, V -> Pt, Mat);
  1105.     V = V -> Pnext;
  1106.     }
  1107.     while (V != NULL && V != VHead);
  1108.  
  1109.     V = VHead;
  1110.     do {
  1111.     Vnext = V -> Pnext;
  1112.     Dx = ABS(V -> Pt[0] - Vnext -> Pt[0]);
  1113.     Dy = ABS(V -> Pt[1] - Vnext -> Pt[1]);
  1114.     Pt[0] = (V -> Pt[0] + Vnext -> Pt[0]) / 2.0;/* Prepare middle point. */
  1115.     Pt[1] = (V -> Pt[1] + Vnext -> Pt[1]) / 2.0;
  1116.     Pt[2] = V -> Pt[2];
  1117.     /* If Dx > Dy fire ray in +Y direction, otherwise in +X direction    */
  1118.     /* and if number of intersections is even (excluding the given point */
  1119.     /* itself) then that direction is the outside, otherwise, its inside.*/
  1120.     Axes = (Dx > Dy ? 1 : 0);
  1121.     if (CGPolygonRayInter(Pl, Pt, Axes) % 2 == 0)
  1122.     {
  1123.         /* The amount we move along Axes is not of a big meaning as long */
  1124.         /* as it is not zero, so MAX(Dx, Dy) guarantee non zero value... */
  1125.         Pt[Axes] -= MAX(Dx, Dy);
  1126.     }
  1127.     else {
  1128.         Pt[Axes] += MAX(Dx, Dy);
  1129.     }
  1130.  
  1131.     /* Now Pt holds point which is in the inside part of vertex V, Vnext.*/
  1132.     /* Put it in the pseudo inside polygon PlIn:                 */
  1133.     if (VInHead) {
  1134.         VIn -> Pnext = AllocVertex(0, 0, NULL, NULL);
  1135.         VIn = VIn -> Pnext;
  1136.     }
  1137.     else {
  1138.         PlIn -> V = VInHead = VIn = AllocVertex(0, 0, NULL, NULL);
  1139.     }
  1140.     PT_COPY(VIn ->Pt, Pt);
  1141.  
  1142.     V = Vnext;
  1143.     }
  1144.     while (V != NULL && V != VHead);
  1145.     VIn -> Pnext = VInHead;
  1146.  
  1147.     MyFree((char *) Pl, ALLOC_POLYGON);/* Free copied (and trans.) vrtx list.*/
  1148.  
  1149.     /* Transform PlIn to the plane where original Pl is... */
  1150.     if (!MatInverseMatrix(Mat, Mat))         /* Find the inverse matrix. */
  1151.     FatalError("GenInsidePoly: Inverse matrix does not exits");
  1152.     VIn = VInHead;
  1153.     do {
  1154.     MatMultVecby4by4(VIn -> Pt, VIn -> Pt, Mat);
  1155.     VIn = VIn -> Pnext;
  1156.     }
  1157.     while (VIn != NULL && VIn != VInHead);
  1158.  
  1159.     return PlIn;
  1160. }
  1161.  
  1162. /*****************************************************************************
  1163. *   Routine to create a polygon out of a list of 4 vertices V1/2/3/4         *
  1164. * The fifth vertex is inside (actually, this is not true, as this point will *
  1165. * be in the positive part of the plane, which only locally in the object...) *
  1166. * the object, so the polygon normal direction can be evaluated uniquely.     *
  1167. *  No test is made to make sure the 4 points are co-planar...             *
  1168. *  The points are placed into a linear line in order.                 *
  1169. *****************************************************************************/
  1170. static PolygonStruct *GenPolygon4Vrtx(VectorType V1, VectorType V2,
  1171.     VectorType V3, VectorType V4, VectorType Vin, PolygonStruct *Pnext)
  1172. {
  1173.     PolygonStruct *PPoly;
  1174.     VertexStruct *V;
  1175.  
  1176.     PPoly = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), Pnext);
  1177.     PT_COPY(V -> Pt, V1);
  1178.  
  1179.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  1180.     PT_COPY(V -> Pt, V2);
  1181.  
  1182.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  1183.     PT_COPY(V -> Pt, V3);
  1184.  
  1185.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  1186.     PT_COPY(V -> Pt, V4);
  1187.  
  1188.     V -> Pnext = PPoly -> V;           /* Make the Vertex list circular. */
  1189.  
  1190.     UpdatePolyPlane(PPoly, Vin);           /* Update plane equation. */
  1191.  
  1192.     SET_CONVEX_POLY(PPoly);               /* Mark it as convex polygon. */
  1193.  
  1194.     return PPoly;
  1195. }
  1196.  
  1197. /*****************************************************************************
  1198. *   Routine to create a polygon out of a list of 3 vertices V1/2/3         *
  1199. * The forth vertex is inside (actually, this is not true, as this point will *
  1200. * be in the positive part of the plane, which only locally in the object...) *
  1201. * the object, so the polygon normal direction can be evaluated uniquely.     *
  1202. *  The points are placed into a linear line in order.                 *
  1203. *****************************************************************************/
  1204. static PolygonStruct *GenPolygon3Vrtx(VectorType V1, VectorType V2,
  1205.             VectorType V3, VectorType Vin, PolygonStruct *Pnext)
  1206. {
  1207.     PolygonStruct *PPoly;
  1208.     VertexStruct *V;
  1209.  
  1210.     PPoly = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), Pnext);
  1211.     PT_COPY(V -> Pt, V1);
  1212.  
  1213.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  1214.     PT_COPY(V -> Pt, V2);
  1215.  
  1216.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  1217.     PT_COPY(V -> Pt, V3);
  1218.  
  1219.     V -> Pnext = PPoly -> V;           /* Make the Vertex list circular. */
  1220.  
  1221.     UpdatePolyPlane(PPoly, Vin);           /* Update plane equation. */
  1222.  
  1223.     SET_CONVEX_POLY(PPoly);               /* Mark it as convex polygon. */
  1224.  
  1225.     return PPoly;
  1226. }
  1227.  
  1228. /*****************************************************************************
  1229. *   Routine to preper transformation martix to do the following (in this     *
  1230. * order): scale by Scale, rotate such that the Z axis is in direction Dir    *
  1231. * and then translate by Trans.                             *
  1232. *    Algorithm: given the Trans vector, it forms the 4th line of Mat. Dir is *
  1233. * used to form the second line (the first 3 lines set the rotation), and     *
  1234. * finally Scale is used to scale first 3 lines/columns to the needed scale:  *
  1235. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord *
  1236. *                |  Bx  By  Bz  0 |  system into T, N & B as required and    *
  1237. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |  then translate it to C. T, N, B are     *
  1238. *                |  Cx  Cy  Cz  1 |  scaled by Scale.                 *
  1239. * N is exactly Dir (unit vec) but we got freedom on T & B - T & B must be on *
  1240. * a plane perpendicular to N and perpendicular between them but thats all!   *
  1241. * T is therefore selected using this (heuristic ?) algorithm:             *
  1242. * Let P be the axis of which the absolute N coefficient is the smallest.     *
  1243. * Let B be (N cross P) and T be (B cross N).                     *
  1244. *****************************************************************************/
  1245. static void GenTransformMatrix(MatrixType Mat, VectorType Trans,
  1246.                         VectorType Dir, RealType Scale)
  1247. {
  1248.     int i, j;
  1249.     RealType R;
  1250.     VectorType DirN, T, B, P;
  1251.     MatrixType TempMat;
  1252.  
  1253.     PT_COPY(DirN, Dir);
  1254.     PT_NORMALIZE(DirN);
  1255.     PT_CLEAR(P);
  1256.     for (i = 1, j = 0, R = ABS(DirN[0]); i < 3; i++) if (R > ABS(DirN[i])) {
  1257.     R = DirN[i];
  1258.     j = i;
  1259.     }
  1260.     P[j] = 1.0;/* Now P is set to the axis with the biggest angle from DirN. */
  1261.  
  1262.     VecCrossProd(B, DirN, P);                  /* Calc the bi-normal. */
  1263.     VecCrossProd(T, B, DirN);                /* Calc the tangent. */
  1264.  
  1265.     MatGenUnitMat(Mat);
  1266.     for (i = 0; i < 3; i++) {
  1267.     Mat[0][i] = T[i];
  1268.     Mat[1][i] = B[i];
  1269.     Mat[2][i] = DirN[i];
  1270.     }
  1271.     MatGenMatScale(Scale, Scale, Scale, TempMat);
  1272.     MatMultTwo4by4(Mat, TempMat, Mat);
  1273.  
  1274.     MatGenMatTrans(Trans[0], Trans[1], Trans[2], TempMat);
  1275.     MatMultTwo4by4(Mat, Mat, TempMat);
  1276. }
  1277.  
  1278. /*****************************************************************************
  1279. *   Routine to update the Plane equation of the given polygon such that the  *
  1280. * Vin vertex will be in the positive side of it.                 *
  1281. *   It is assumed the polygon has at list 3 points...                 *
  1282. *****************************************************************************/
  1283. void UpdatePolyPlane(PolygonStruct *PPoly, VectorType Vin)
  1284. {
  1285.     int i;
  1286.     VectorType V1, V2;
  1287.     RealType Len;
  1288.     VertexStruct *V;
  1289.  
  1290.     V = PPoly -> V;    PT_SUB(V1, V -> Pt, V -> Pnext -> Pt);
  1291.     V = V -> Pnext;    PT_SUB(V2, V -> Pt, V -> Pnext -> Pt);
  1292.  
  1293.     VecCrossProd(PPoly -> Plane, V1, V2);           /* Find Plane Normal. */
  1294.     /* Normalize the plane such that the normal has length of 1: */
  1295.     Len = PT_LENGTH(PPoly -> Plane);
  1296.     for (i = 0; i < 3; i++) PPoly -> Plane[i] /= Len;
  1297.  
  1298.     PPoly -> Plane[3] = (-DOT_PROD(PPoly -> Plane, PPoly -> V -> Pt));
  1299.  
  1300.     if (DOT_PROD(PPoly -> Plane, Vin) + PPoly -> Plane[3] < 0) {
  1301.     /* Flip plane normal and reverse the vertex list. */
  1302.     ReverseVrtxList(PPoly);
  1303.     for (i = 0; i < 4; i++) PPoly -> Plane[i] = (-PPoly -> Plane[i]);
  1304.     }
  1305. }
  1306.  
  1307. /*****************************************************************************
  1308. *   Routine to update the Vertex Pt normal equation. The normal should point *
  1309. * InPt but should be perpendicular to the line Pt-PerpPt if Perpendicular is *
  1310. * TRUE. THe normal is normalized to a unit length.                 *
  1311. *****************************************************************************/
  1312. static void UpdateVertexNormal(NormalType Normal, PointType Pt, PointType InPt,
  1313.                     int Perpendicular, PointType PerpPt)
  1314. {
  1315.     VectorType V1, V2;
  1316.  
  1317.     PT_SUB(Normal, InPt, Pt);
  1318.  
  1319.     if (Perpendicular) {
  1320.     PT_SUB(V1, PerpPt, Pt);
  1321.     VecCrossProd(V2, V1, Normal);
  1322.     VecCrossProd(Normal, V2, V1);
  1323.     }
  1324.  
  1325.     PT_NORMALIZE(Normal);
  1326. }
  1327.